Welcome

The main focus is on two areas - experimental design under treatment interference and multi-objective experimental design - that are active research areas.

Common DoE assumptions

To a greater or lesser degree, the ideas and principles from Weeks 1-3 are dependent on a number of assumptions holding:

• that the expected response from a given unit only depends on the treatment applied to that unit, and not on the treatments applied to any other units;

• for factorial, response surface and optimal designs, that a reasonable approximating statistical model can be specified;

• for optimal designs, that the aim of the experiment can be neatly encapsulated in a single mathematical expression or objective function.

Relaxing the assumptions

This week, we will focus on approaches that allow us to relax one or the other of these assumptions:

• in Part 1 (Experiments with interference), we will introduce methods for designing and analysing experiments when treatment interference is anticipatedl

• in Part 2 (Multi-objective experimentation), we will introduce multi-objective (compound) design optimality criteria that address multiple experimental aims simultaneously.

SUTVA

Standard model for analysing a designed experiment,

y_{i} = \mu + \tau_{r(i)} + \varepsilon_i\,, \tag{1}

with the aim of estimating treatment differences \tau_j - \tau_k. Here, r(i) \in \{1,\ldots,t\} indicates which treatment was allocated to the ith unit (i = 1,\ldots,n).

This model makes the stable unit treatment value assumption (SUTVA), which states that the response from any particular unit is unaffected by the assignment of treatments to other units (Cox 1958, sec. 2.4).

Ignoring substantial treatment interference, as in model (1), can lead to biased estimates of differences between the direct treatment effects t_r.

Mitigate of treatment interference may be possible, e.g., by adding “wash-out” periods in the clinical trial or “guard plots” in the agricultural experiment.

But in many cases this may not be possible (or ethical) or there may be interest in the indirect effect of each treatment; for example, the viral effect of the adverts in the marketing experiment.

Hence, it is of interest to study designs and models which account for treatment interference.

Cross-over trials

In a cross-over trial, each subject is assigned a sequence of treatments across different time periods. Interest is in comparing individual treatments, not sequences, and the experimental units are the periods within each subject.

Cross-over trials are common in studies of chronic conditions, where repeated treatment is required.

Advantages include

• each subject acts as their own control;
• and treatment comparisons can be made “within subject”.

See Jones and Kenwood (2015) and Bose and Dey (2015).

Order and carry-over effects

However, this feature of within subject comparison can also bring disadvantages. Principally,

• the order in which the treatments are applied may impact the outcome;
• there may be treatment interference between periods of the design (typically called carry-over in a cross-over trial). That is, the outcome from period i may depend on the treatment applied in period i-1 (first-order interference) or even from earlier periods.

2x2 cross-over trial

The simplest form of cross-over design concerns t=2 treatments and p=2 periods.

There are clearly only two possible sequences and each subject in the trial is randomised to one of the two. A wash-out period is may be inserted between the two treatment periods.

Example

Investigate the efficacy of an inhaled drug (A), compared to a control (B), for patients suffering from chronic obstructive pulmonary disease (COPD). The response was the mean expiratory flow rate (PEFR) based on readings recorded each morning by the subjects.

Linear model for cross-over trials

The traditional linear model used with cross-over experiments contains terms corresponding to subjects, period, treatment and interference:

\begin{split} y_{ij} = \mu + \alpha_i + \beta_j + \tau_{r(i,j)} + \rho_{r(i-1,j)} + \varepsilon_{ij}\,,\\ i = 1,\ldots,p;\, j = 1,\ldots, n\,, \end{split} \tag{2}

• r(i,j) \in \{1,\ldots,t\} denotes the treatment allocated to the jth subject in the ith period,
• \alpha_i is the ith period effect,
• \beta_j is the jth subject effect,
• \rho_{r(i-1,j)} is the interference (indirect, carry-over) effect, with \rho_{r(0,j)}=0.

Analysis of variance

Notes:

• the between sequence sums of squares tests \rho_1 = \rho_2;
• the correct demominator sums of squares is between subjects;
• we need the assumption of no interference to test for direct treatment effects (as we cannot have complete balance as a treatment does not follow itself).

Larger cross-over designs

A balanced design has each treatment occuring the same number of times in each period, and each treatment following every other treatment the same number of times, with no treatment following itself.

• A 2\times 2 trial is balanced.
• Every pair of estimated direct treatment differences, \hat{\tau}_k = \hat{\tau}_l has the same variance.

A strongly balanced design has every treatment followed by every other treatment, including itself.

• Direct and indirect effects are estimated independently, simplifying the analysis and interpretation, and lowering the variance of estimators of the indirect effects.

Balanced 4x4 Latin square

Strongly balance extra-period design

Networked experiments

Experiments from outside the clinical arena can also violate SUTVA.

For example, online controlled experiments on websites and social media platforms (Larsen et al. 2024).

• Connections between users can lead to treatment interference.
• There could be distinct communities within the network, with more similar responses expected from users in the same community.

Figure 1: Subset of the Facebook network from the Stanford Data Collection, Colours inducated 24 distinct blocks, or communities, of users.

Other examples

Graphs of designs

• Vertices and edges
• Adjacency matrix
• Blocks

Suppose that n units are formed into a network with connections representing possible treatment interference.

This network can be represented as a graph \mathcal{G} = (\mathcal{V}, \mathcal{E}),

• vertex set \mathcal{V} represents units;
• edge set \mathcal{E}, of size l, represents the connections.

Connections can be succinctly represented via the adjacency matrix.

• A = [A]_{jh}, an n\times n matrix with A_{jh} = \in [0, 1]
• Undirected graphs have A_{jh} = A_{hj}.
• Lack of an edge between nodes j and h in \mathcal{E} leads to A_{jh} = A_{hj} = 0.

	A	B	C	D	E
A	0	1	0	1	1
B	1	0	1	0	0
C	0	1	0	1	0
D	1	0	1	0	1
E	1	0	0	1	0

For some applications, necessary blocking factors may be obvious or based on covariates external to the graph, e.g., age or sex.

For others, it may be necessary or desirable to base the blocks on the graph structure itself, e.g., using spectral clustering (Koutra, Gilmour, and Parker 2021).

Linear network model

The adjacency matrix can be used to incorporate indirect treatment effects into a model for the experiment Koutra, Gilmour, and Parker (2021):

\begin{split} y_{ij} = \mu + \beta_i + \tau_{r(i,j)} + \sum_{g=1}^{b}\sum_{h=1}^{n_{g}} A_{\left\{ij,gh\right\}}\gamma_{r(g,h)} +\varepsilon_{ij}\,, \\ \quad i=1,\ldots,b\,,\,j = 1,\ldots, n_i\,. \end{split} \tag{3}

• \beta_i is the ith block effect;
• \tau_k and \gamma_k are the direct and indirect treatment effects

Example

For the graph in Figure Figure 2, with adjacency matrix in Table 7, the response from node A, the first unit in block 1, would be modelled as:

y_{11} = \mu + \beta_1 + \tau_{r(1,1)} + \gamma_{r(1,2)} + \gamma_{r(1,3)} + \gamma_{r(2,1)} + \varepsilon_{11}\,,

with the indirect treatment effects resulting from the edges between node A and nodes D and E (in block 1) and node B (in block 2). The linear network effects model can be estimated using least squares or maximum likelihood.

A-optimality

For efficient estimation of direct treatment differences, designs are sought that minimise the average variance of the pairwise differences:

\phi_{\tau}=\frac{2}{t(t-1)} \sum_{s=1}^{t-1}\sum_{s'=s+1}^t \text{var}(\widehat{\tau_s-\tau_{s'}})\,. \tag{4}

Similarly, we can define a criterion for efficient estimation of indirect treatment differences, that minimises

\phi_{\gamma}=\frac{2}{t(t-1)} \sum_{s=1}^{t-1}\sum_{s'=s+1}^t \text{var}(\widehat{\gamma_s-\gamma_{s'}})\,. \tag{5}

Designs can be found via application of standard optimisation algorithms, such as point exchange (Cook and Nachtsheim 1980).

Example - co-authorship network

Links between academics within a university research group (Koutra, Gilmour, and Parker 2021).

• 22 nodes, split into three blocks, and 27 edges.
• Interest lies in estimating the effects of two treatments.

Figure 3: Block designs for a co-authorship network with colours indicating blocks (identified via spectral clustering) and plotting symbol indicating allocation to treatment 1 or 2. Left: optimal design for estimation of direct effects. Right: optimal design for estimation of indirect effects.

Notes

Estimating direct effects

• The \phi_{\tau}-optimal design is balanced, with equal replication of each treatment.
• Treatment allocation is also balanced within each block.
• Nodes allocated to each treatment have similar first- and second-order degrees

Estimating indirect effects

• Treatment allocation in the \phi_{\gamma}-optimal design is highly dependent on the network.
• Highly connected nodes tend to receive a different treatment from their surrounding, less connected, nodes.

Comparisons to other designs

Quantitative comparisons can be made between the block network designs (BNDs) from Figure 3 and the optimal designs that would result from models that

• ignore blocks and network structure (CRD: completely randomised design);
• ignore network structure (RBD: randomised block design);
• ignore blocks (LND: linear network design).

Efficiencies for estimating direct effects

Efficiencies are calculated within row.

Efficiencies for estimating indirect effects

The loss of efficiency for designs that ignore network structure is now large (>80%).

Case study - agricultural experiment

The impact of neighbouring plots in field trials has been widely considered, including through study of indirect treatment effects (Besag and Kempton 1986).

Figure 4: Example field layout, as used at Rothamsted

A typical layout of a field trial is shown in Figure 4, clearly showing the proximity of neighbouring plots.

Row-column experiment

Experiments conducted at Rothamsted to study the differences in natural cereal aphid colonization (Koutra et al. 2023).

• 21 different wheat varieties.
• Units arranged in a 14\times 6 grid of 1m x 1m plots.
• There are sufficient plots for each treatment to be replicated four times.
• Data from 2016 experiment.

Layout and 2016 experiment

Figure 5: Plot layout for the agricultural example with treatment allocation from the 2016 design. Numbers indicated treatments allocated to each plot.

Treatment interference

Treatment interference was thought possible due to the differing levels of susceptibility of different varieties and the strong possibility of aphids moving from plot to plot.

Differing structures governing this interference were considered, represented as graphs.

Figure 6: Network and optimal design for the wheat field trial. Numbers indicated treatments allocated to each plot.

Statistical model

In addition to direct and indirect treatment effects, the analysis of the experiment needed to account for the spatial structure through the inclusion of blocking factors and row-column effects.

\begin{aligned} y_j &=\mu+\tau_{r(j)}+ R_i+C_k+(RC)_{ik}+r_{ig}+c_{kh} \\ & +\left(rC\right)_{igk}+\left(Rc\right)_{ikh} +\sum_{j'} A_{jj'} \gamma_{r(j')} + \varepsilon_{j}\,, \end{aligned} \tag{6}

with R, C and RC representing the effects of super-rows and super-columns, and their interaction (super-blocks). Effects r and c are of rows and columns nested inside super-blocks.

Analysis of variance from 2016 experiment

Optimal design for estimating direct effects

Figure 7: Network and optimal design for the wheat field trial. Numbers indicated treatments allocated to each plot.

Notes and comparisons

• There is a good spatial spread of treatments;
• but it is also quite common for pairs of connected units to share a treatment.

Both these features have been observed in previous row-column and network designs; see Freeman (1979), Parker, Gilmour, and Schormans (2016) and Koutra, Gilmour, and Parker (2021).

A comparison to designs found under different models shows a big loss in efficiency. The efficiency of the 2016 design (a resolvable row-column design) was 0.5.

Response surface designs

We will discuss such methods in the context of response surface models of the form

\begin{split} Y_i & = \beta_0 + \sum_{j = 1} ^ {p} \beta_j x_{ij} + \varepsilon_{i} \\ & = \beta_0 + \boldsymbol{x}_{i1}^{\mathrm{T}}\boldsymbol{\beta}_1 + \varepsilon_{i}\,, \end{split} \tag{7}

• \beta_0 and \boldsymbol{\beta}_1^{\mathrm{T}} = (\beta_1, \ldots, \beta_{p}) are unknown parameters.
• \boldsymbol{x}_{i1}^{\mathrm{T}} = (x_{i1}, \ldots, x_{ip}) holds values of the p predictors for the ith run.

Standard optimality criteria

The most common design selection criteria aim at estimation for model (7) with \sigma^2 assumed known.

D-optimality \phi_{D}(\mathcal{D}) = \left|\left[X_1^{\mathrm{T}}\left(I_n - \frac{1}{n}J_n\right)X_1\right]^{-1}\right|\,.\\ L-optimality \phi_L(\mathcal{D}) = \text{tr}\left\{L^{\mathrm{T}}\left(X_1^{\mathrm{T}}(I_n - \frac{1}{n}J_n)X_1\right)^{-1}L\right\}\,.

Here, we treat the intercept \beta_0 as a nuisance parameter.

Both criteria assume model (7) is correctly specified.

Primary and potential terms

Box and Draper (1959) introduced the idea of discrepancy between the assumed response surface model and an encompassing “true” model:

\begin{split} Y_i & = \beta_0 + \sum_{j = 1} ^ {p} \beta_j x_{ij} + \sum_{j = p+1} ^ {p+q} \beta_j x_{ij} +\varepsilon_i \\ & = \beta_0 + \boldsymbol{x}_{i1}^{\mathrm{T}}\boldsymbol{\beta}_1 + \boldsymbol{x}_{i2}^{\mathrm{T}}\boldsymbol{\beta}_2 + \varepsilon_i\,, \end{split} \tag{8}

where \boldsymbol{x}_{i2}^{\mathrm{T}} = (x_{i(p+1)}, \ldots, x_{i(p+q)}) holds the additional q polynomial terms, with associated parameters \boldsymbol{\beta}_2^{\mathrm{T}} = (\beta_{p+1}, \ldots, \beta_{p+q}).

DuMouchel and Jones (1994) labelled the polynomial terms in the assumed model as primary and the additional terms in the encompassing model as potential.

Mean squared error

One desirable aim is to be able to estimate \boldsymbol{\beta}_1 from model (7) protected from contamination from the potential terms.

Define the MSE matrix for \hat{\boldsymbol{\beta}} (Montepiedra and Fedorov 1997):

\begin{split} \text{MSE}\left(\hat{\boldsymbol{\beta}}_1\right)& = \mathtt{E}_{\boldsymbol{Y}}[(\hat{\boldsymbol{\beta}}_1 -\boldsymbol{\beta}_1)(\hat{\boldsymbol{\beta}}_1 - \boldsymbol{\beta})_1^\top]\\ & = \sigma^2[X_1^{\mathrm{T}} (I_n - \frac{1}{n}J_n) X_1]^{-1} + A_1\boldsymbol{\beta}_2\boldsymbol{\beta}_2^{\mathrm{T}} A_1^{\mathrm{T}}\,, \end{split} \tag{9}

where

A_1 = \left[X_1^{\mathrm{T}} \left(I_n - \frac{1}{n}J_n\right) X_1\right]^{-1}X_1^{\mathrm{T}} \left(I_n - \frac{1}{n}J_n\right)X_2

is the p\times q alias matrix between the primary and potential terms (excluding the intercept).

MSE(L) optimality

An analogy of variance-based alphabetic criteria is to consider functionals of this matrix.

\begin{split} \phi_{MSE(L)}(\mathcal{D}) & = E\left\{\text{trace}\left[\text{MSE}\left(\hat{\boldsymbol{\beta}}_1\right)\right]\right\} \\ & = \text{trace}\left\{E\left[\text{MSE}\left(\hat{\boldsymbol{\beta}}_1\right)\right]\right\} \\ & = \text{trace}\left[\sigma^2M^{-1} + E\left(A_1\boldsymbol{\beta}_2\boldsymbol{\beta}_2^\top A_1^\top\right)\right] \\ & = \sigma^2\text{trace}\left[M^{-1} + \tau^2 A_1^\top A_1\right]\,. \end{split} The expectation is taken with respect to a normal prior distribution for \boldsymbol{\beta}_2\sim \mathcal{N}\left(\boldsymbol{0}_q, \sigma^2\tau^2I_q\right) for \tau^2>0.

Pure error criteria

When uncertainty about the assumed model is being acknowledged, it is important that sufficient pure error degrees of freedom exist in the design to provide an unbiased estimator for \sigma^2.

Gilmour and Trinca (2012) suggested a class of criteria that explicitly incorporate the F-distribution quantiles on which parameter confidence regions depend.

DP-optimality \phi_{(DP)_S}(\mathcal{D}) = F_{p,d;1-\alpha}^p\phi_{D}(\mathcal{D})\,.

LP-optimality \phi_{LP}(\mathcal{D}) = F_{1,d;1-\alpha}\phi_L(\mathcal{D})\,.

where d = n-t is the number of replicated treatments in the experiment, \alpha is a pre-chosen significance level and F_{df1, df2; 1-\alpha} is the quantile of an F-distribution with df1 and df2 degrees of freedom such that the probability of being less than or equal to this quantile is 1-\alpha.

Model sensitivity

The ability of the design to make inference about the potential terms, and hence detect any lack of fit in the direction of model (8) can be quantified via functionals of R + \frac{1}{\tau^2}I_q, which is proportional to the posterior variance for \boldsymbol{\beta}_2.

LoF-DP-optimality \phi_{LoF-DP}(\mathcal{D}) = F^q_{q, d; 1-\alpha_{L}} \left|R + \frac{1}{\tau^2}I_q\right|^{-1}\,. LoF-LP-optimality \phi_{LoF-LP}(\mathcal{D}) = F_{1, d; 1-\alpha_{L}} \text{tr}\left\{L^\top\left(R + \frac{1}{\tau^2}I_q\right)^{-1}L\right\}\,.

Both criteria target designs with matrices X_1 and X_2 being (near) orthogonal to each other, which will also maximise the power of the lack-of-fit test for the potential terms.

MOODE

Multi-objective optimal design of experiments can be achieved via a compound criterion objective function constructed via a weighted product of individual objective functions.

Egorova and Gilmour (2023) defined a trace-based compound criteria as

\phi_{trace}(\mathcal{D}) = \phi_{LP}(\mathcal{D})^{\kappa_{LP}}\times \phi_{LoF-LP}(\mathcal{D})^{\kappa_{LoF-LP}} \times \phi_{MSE(L)}(\mathcal{D})^{\kappa_{MSE(L)}}\,,

with all weights \kappa \ge 0 and \kappa_{LP} + \kappa_{LoF-LP} + \kappa_{MSE(L)} = 1.

These compound criteria, along with their componenet criteria, are implemented in the R package MOODE (Koutra et al. 2024), available on CRAN.

Example

The 12-run Plackett-Burman design is perhaps the most widely used, and studied, non-regular fractional factorial design.

• Orthogonal estimation of the main effects of up to 11 two-level factors
• Main effect estimator for each factor is biased by all two-factor interactions not involving that factor.

We find alternative two-level designs for k=3,\ldots,9 factors using the trace-based compound criterion under five different sets of criteria weights.

The primary model consists of all k main effects, with the potential model also including all two-factor interactions.

The MOODE package can be used to find designs under these models and criteria.

Weights

Results

Figure 8: Efficiencies for MOODE and other optimal designs under different criteria.

Results

Notes 1

LP-efficiency:

• L-optimal designs lack replication, few pure error degrees of freedom and zero efficiency under LP-criterion.
• Similar for MSE(L)-optimal designs for larger k.
• Compound optimal designs at least 50% LP-efficient, usually more.

MSE(L)-efficiency:

• LP-optimal and compound designs have too much replication, and low efficiency for small k.
• For larger k, less replication is possible and efficiency improves

Notes 2

For k=4:

• LP-optimal design has only 5 distinct points (the minimum number possible),
• Compound desing has 8 distinct points
• MSE(L)-optimal design has 12 distinct points (obviously the maximum)

None of these designs are orthogonal in the main effects, a property of the Plackett-Burman design that is compromised to obtain either

• PE degrees of freedom (LP);
• or better robustness from the potential terms (compound and MSE(L)).

Compound and MSE(L)-optimal designs achieve orthogonality between the main effects and two-factor interactions, i.e. A_1 is a zero matrix.

Conclusions

Two areas of active research:

• treatment inference;
• multi-objective design.

Both reduce the reliance on assumptions that may be unrealistic in many cases.

The topics could also be combined. e.g., multi-objective designs could be sought for networked experiments to estimate both direct and indirect treatment effects.

Further work

Both topics also intersect with other research areas in design of experiments.

• Networked experiments is also an active area within the causal inference community (e.g., Hudgens and Halloran 2008); a workshop was held at King’s in the summer of 2024.

• Increasingly, experiments are taking place on very large networks, particularly online experimentation e.g., on social media (e.g., Nandy et al. 2020); connections can be made to methods for subsampling large data using design of experiments principles, e.g., Yu, Ai, and Ye (2024).